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\ ABSTRACT 

\ 

Based  on  the  assumption  that  the  wave  energy  associated  with  a narrow 

frequency  band  stayed  within  the  band  upon  refraction,  a numerical  model  based 
on  finite  difference  method  has  been  developed.  This  model,  utilizing  electronic 
computer,  proves  to  be  quite  convenient  to  determine  shallow  water  wave  spectra 
at  designated  spatial  locations  of  irregular  bottom  contours  provided  the 
deepwater  wave  spectrum  ic  given.  The  model  has  been  verified  against  analytical 
solution  for  case  of  two-di.nensional  parallel  bottom  contours  and  compared  in  good 
agreement  with  field  measurement  at  location  well  beyond  the  surf  zone. 
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CHAPTER  I INTRODUCTION 


The  needs  for  more  realistic  representation  and  accurate  prediction  of 
sahllow  water  environment  are  increasing  for  both  scientific  and  engineering 
purposes,  in  particular,  due  to  the  advent  of  anticipated  continental  shelf 
development.  It  is  also  an  accepted  opinion  that  wave  spectrum  which  utilizes  the 
ensemble  of  random  surface  oscillations  provides  a more  realistic,  and  perhaps  more 
accurate,  representation  of  ocean  waves  than  monochromatic  wave  train  of  regular 
oscillation.  In  the  past,  considerable  efforts  have  been  devoted  in  the  development 
of  deep  water  wave  spectrum  to  the  extent  that  one  has  reasonable  confidence  on 
its  applicability  to  describe  wind  generated  ocean  waves.  The  goal  of  this  work 
Is  to  develop  a practical  method  of  transforming  this  deepwater  water  wave  spectra 
into  shallow  water  of  irregular  bottom  topography. 


A number  of  investigators  (I.onguet-Higgins , 1956;  Collins,  1970;  etc.) 
have  studied  the  theoretical  formulation  of  the  transformation  of  wind  wave 
spectrum  in  shoaling  water  by  considering  the  conservation  of  energy  between 
adjacent  wave  rays  much  the  same  way  of  the  wave  refraction  in  monocromntic  wave 
train.  Based  on  Longuet-Higgins'  formulation,  Knrlesson  (1969)  developed  a 
method  to  compute  spectrum  transformation  over  parallel  bottom  contours.  Collins 
(1972)  extended  the  work  to  the  inclusion  of  bottom  friction  and  have  mentioned 
the  application  to  the  irregular  bottom  topography.  However,  his  numerical  scheme 
traces  wave  energy  along  wave  rays  which  makes  the  computational  procedure  quite 
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impractical  in  determining  shallow  water  wave  spectrum  at  designated  locations 
over  an  area  of  irregular  offshore  bottom  configuration.  More  recently, 

Krasitkiy  (1974)  has  derived  explicit  analytical  solutions  for  the  spectrum  trans- 
formation over  two-dimensional  parallel  bottom  contrours. 

The  present  study  presents  a numerical  method  which  proves  to  be  effective 
to  transform  deepwater  wave  spectrum  into  shallow’  water  regions.  It  is  based 
on  the  assumption  that  the  wave  energy  associated  with  a narrow  frequency  band 
stays  within  this  band  upon  refraction  when  waves  propagate  ftom  deep  into 
shallow  waters.  A recently  developed  wave  refraction  program  (Noda,  et_  aJL . , 1974) 
is  modified  to  compute  Che  spatial  modification  of  the  incident  wave  energy 
within  the  designated  frequency.  The  shallow  water  wave  spectrum  at  any  specified 
location,  is  then  simply  the  contributions  resulting  from  energy  associated 

I 

with  different  frequency  bands  of  the  spectrum. 

Examples  of  the  spectrum  distribution  over  parallel  bottom  contours,  irregular 
bottom  bathymetry  and  rhythmic  topography  were  illustrated.  Comparisons  of  the 


CHAPTER  II  TRANSFORMATION  OF  ENERGY  SPECTRUM 


Let  F(k,r,t)  be  the  spectral  density  function  defined  as  energy  per  unit 
area  per  unit  wave  number  space,  the  total  energy  is  equal  to  the  summation  of 
the  spectral  density  over  the  whole  wave  number  space,  i.e., 

(1) 


E = / F(k,r, t)dk 

J 0 


where  k - magnitude  of  the  wave  number  = 2,;r/L 

L = wave  length 
-> 

r = position  vector 
t = time 


For  the  case  of  steady  state,  the  wave  energy  contained  within  a small  increment 
dk  can  be  expressed  as 
' k+dk 
k 

Through  Taylor's  expansion  and  neglecting  higher  orders,  we  have 


AE 


f k+dk 

= j F(C,?)dC 


AE  = F(k,r)dk 


(2) 


It  is  known  that  for  linear  wave  svstems,  the  following  dispersion  relationship 
exists 

2 

w = gk  Tanh  kh  (3) 

where  to  = wave  angular  frequency  = 2tt/T 
T = wave  period 
h = local  water  depth 
g = gravitational  acceleration 
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Differentiating  the  above  equation  results  in 
dw  = Cg  dk 

where  Cg  = wave  group  velocity  = nc 
n * (1  + 2kh/sinh  2kh)/2 

c = wave  phase  velocity  = L/T 


(A) 


Utilizing  the  above  relation,  the  wave  number  spectral  density  function  F(k,r) 
can  be  transformed  into  polar  frequency  direction  space  as 


F(k,r)  = jp-  4>(u,  a,  r) 


where  <f>  = the  frequency  directional  wave  spectral  density  function 
r = magnitude  of  space  vector  "r 
a = polar  angle  of  the  spectral  vector 
Combining  Equations  (2)  and  (5),  one  obtains 

k^E 


(5) 


4>(w,  a,r) 


d w 


(6) 


Assuming  that  the  wave  energy  associated  with  a narrow  frequency  band  would  stay 
within  this  band  upon  wave  refraction,  a relationship  between  shallow  water  and 
deepwater  wave  energy  spectra  can  be  obtained  by 


9 k A£  o 
o 0 

where  the  subscript  o denotes  deepwater  conditions. 

For  small  amplitude  waves,  it  is  known  that  (Lamb,  1932) 

where  Y = specific  weight  of  the  fluid 
H = wave  height 


(7) 


It  is  plausible  to  assume  here  that 


(8) 


where  AH  can  be  interpreted  as  the  contribution  to  the  total  wave  height  due  to 
the  associated  increment  of  wave  energy.  Based  on  this  relationship,  Equation  (7) 
can  be  rewritten  as 


v k . .2  ^o 

o (AHq) 


(9) 


This  equation  which  clearly  reveals  the  effects  of  shoaling  and  dispersion 
is  the  key  equation  for  wave  energy  transformation  computations.  It  is  evident 
from  this  equation  that  the  task  of  computing  <j>  requires  the  computation  of 
wave  height  and  wave  number  transformations  at  designated  locations. 


CHAPTER  III  TRANSFORMATION  OF  WAVE  CLIMATE 


To  compute  wave  height  and  wave  number  transformation  at  designated  locations 
in  shallow  water,  a numerical  program  on  wave  refraction  and  shoaling  recently 
developed  by  Noda  et.al.  (1974)  is  adopted  and  modified.  This  method  like  most 
of  the  existing  techniques  in  computing  wave  refraction,  neglected  the  nonlinear 
effects,  bottom  friction,  turbulent  dissipation  and  bottom  percolation.  The 
procedure  involves  solving  a pair  of  simultaneous-  equations  for  the  wave  direc- 
tion and  the  wave  height.  The  first  equation  is  derived  from  the  fact  that  the 
wave  number  vector  field,  k,  is  irrotational  (Phillips,  1966).  In  Cartesian 
coordinates  as  shown  in  Figure  1,  the  following  equation  results 


„ 39  , . . 30  1 . .3k  . . 3k, 

3x  3y  k 3y  3x 

where  9 = angle  of  wave  incidence  to  the  beach  normal 

x = Cartesian  coordinate  perpendicular  to  the  shoreline 

y = Cartesian  coordinate  parallel  to  the  shoreline 


(10) 


The  partial  derivatives  of  k were  determined  from  the  wave  dispersion  relationship 
expressed  in  Equation  (3).  The  second  equation  is  obtained  from  the  condition 
of  the  steady  state  conservation  of  energy  (Phillips,  1966). 


3C 


(cg  C039)  I If  + (CB  sin6>  I U + cos0  ^ - C sin0 
S n jx  g H ay  3x  g 3x 


3C 


(ID 


+ sin 9 + C cos9  — = 0 
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Ik  the  surf  zone,  waves  break  due  to  hydrodynamic  instability.  In  this 
region  Equation  11  is  no  longer  valid  and  the  wave  height  is  obtained  through 
the  use  of  an  empirical  equation  developed  by  LeMehaute  and  Koh  (1S67). 

- 0.12  Cp)  tanh  (kJO  (12) 

"b 

where  the  subscript,  b,  denotes  breaking  conditions. 

These  sets  of  equations.  Equations  10,  11  and  12,  when  solved  numerically, 
enables  one  to  determine  the  local  wave  height,  wave  number  and  wave  angle  at 
any  designated  location  in  the  area  of  interest  such  as  the  grid  points  in 
Figure  2. 

With  the  local  wave  climate  information  in  hand,  the  wave  energy  spectral 

density  function  can  be  easily  obtained  by  virtue  of  Equation  (9).  It  is 

observed  here  that  the  wave  shoaling  effect  is  purely  a function  of  wave  frequency 

and  bottom  topography.  Therefore,  in  computing  the  ratio  of  AK/AHq  as  required 

in  Equation  (9),  it  is  nonessential  to  know  the  absolute  magnitude  of  AH  provided 

o 

the  cave  does  not  break.  This  fact  considerably  simplifies  the  computation  of 
wave  spectrum  outside  the  surf  zone. 


CHAPTER  IV  NUMERICAL  PROCEDURES 


Computational  Scheme.  The  computation  of  the  transformation  of  wave  climate 

follows  closely  the  method  developed  by  Noda  (1974).  Equations  of  (10)  and  (11) 

are  written  in  finite  difference  form  using  a forward  difference  in  x,  and  a 

backward  difference  in  y for  the  local  grid  system  shown  in  Figure  3.  Each  grid 

was  assigned  a mean  depth  h and  an  area  Ax  Ay.  All  major  quantities — H,  G, 

k — are  computed  at  the  grid  center.  Values  of  sin  0.  . and  cos  0 , are  deter- 

i.J  i*J 

mined  using  an  average  of  these  quantities  from  four  surrounding  grid  blocks 
approximated  to  second  order  in  a Taylor  series  expansion,  i.e., 

sinQ  , , = 7-  (sinQ  ...  . + sin9  . . . + sin0  , ...  + sin0  ... ) 
i,j  4 i+l,j  i-l.J  i.!+l  i.J-1 


+ 8 Uei+l,j  - 8l-l,jKc0sei-l,j  ' cosei+l.j>  + C6l,j+1  - 


(“s6  l,j+l  - “s3  (13) 

etc. 

The  numerical  iterations  are  carried  to  an  acceptable  error,  in  this  case, 

■ - fold 

vJ 

where  X represents  all  the  variables  under  consideration. 

The  computation  of  spectrum  transformation  does  not  involve  numerical 
integration  and  is  straight  forward  from  Equation  (9). 


0.1% 


(14) 
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Figure  2 Grid  System  Used  in  Computation 


Figure  3 Local  Grid  Description 


Boundary  Conditions.  The  condition  at  the  shoreward  boundary  is  simply  that  the 
wave  angles  are  normal  to  the  shoreline  at  the  boundary.  Two  different  options 
are  built  in  the  program  for  the  side  boundaries. 

A.  For  rhythmic  topography — In  this  case  the  wave  condition  repeats  itself 
at  space  intervals  equal  to  the  inceger  multiplications  of  the  periodic  length. 
The  boundary  condition  can  be  self-generated  by  adding  additional  grid  column 
such  that 


*i-  w-i  ■ h.2 

where  N+l  = periodic  length. 

A.  For  slowly  varying  irregular  topography — The  conditions  at  the  upwave  side 
are  computed  by  extending  the  bottom  profiles  two  dimensionally  such  that 
the  wave  refraction  follow  the  Snell  law  up  to  the  boundary.  The  downward  side 
boundary  condition  needs  not  be  prescribed  for  this  case. 

Input  Wave  Condition.  Strictly  speaking,  directional  spectrum  of  the  form 
$ (m,0)  should  be  prescribed  at  the  seaward  boundary.  Since  our  present  know- 
ledge  on  directional  spectrum  is  quite  inadequate,  the  program,  therefore,  is 
written  to  accept  the  deepwater  wave  spectrum  of  the  following  form  as  input: 


$ (w  ,0  ) = S (u>)  K(9  ) 

o 


where  K(0)  is  an  artificial  weighting  function  to  estimate  the  effect  of  direc- 
tional spreading.  For  the  sample  computations  to  be  illustrated  S(w)  is  taken 
to  be  the  Pierson-Moskowitz ' s deep  water  wave  spectrum: 

2 

S(cj)  = a ^ exp  [-3  (-^)^]  (17) 


where  U is  the  wind  speed,  a = 8.1  x 10  and  6 = 0.74. 
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1 


K(0) 


{ 


(13) 


The  function  of  K(0)  assumes  the  form  of 

(8/3*)  cos4  9 | © 1 < */2 

0 |e  | 1 tt/2 

to  facilitate  comparisons  with  analytical  results  of  Krasitskiy  (1974). 

As  the  seaward  boundary  is  often  set  at  finite  water  depths  to  avoid  excessive 
computational  time  the  deep  water  wave  spectrum  is  transformed  to  the  finite 
depths  at  the  boundary  using  analytical  solution  of  Krasitskiy  (1974)  with  the 
assumption  that  the  offshore  topography  can  be  approximated  by  two-dimensional 
parallel  contours. 


■ 
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CHAPTER  V VERIFICATION  AND  APPLICATIONS 


Sample  computations  of  thr.  -a  different  of  bottom  hydrographs  are 

illustrated  her.  t . p rallel  contour 

case  such  that  the  results  can  be  comp;  i . ..  . n of  Krasitskiy 

(1974).  The  calculations  trried  out  for  a wind  speed  of  20.5m/sec  (the 

same  value  chosen  by  krasitskiy)  ar  v : yc  m - • '05.  Figure  4 com- 
pares the  numerical  result  with  the  analytical 

ones  for  water  depths  of  , 25M  and  15M  es  nr  numerical  scheme  yields 

satisfactory  results  equ  . . ■ . :.  r.  from 

the  analytii  il  values  are  not  Inevitable  due 

to  the  finite  grid  size.  . . ..  LI  e th  numerical  errors 

at  the  expense  of  increasing  computational  th  . 

The  second  c ■ is  a i ...  Ison  bi  ui  with  the  field 

measurements.  The  spec  ...  the  Lav  re  coast,  facing 

the  Atlantic  Ocean,  locat  ’d  approxi::  .tc  . y ■ 1 nilt  north  * f Indian  River  Inlet 
(Figure  5)  wa  ■ exan  Lned.  easut  1 t by  pacitanci  -• ype 

wave  gage  mounted  on  i wave  tower  i , ;r  . ; . 4t:.  deep.  The 

location  of  the  tower  and  the  bo,.  a f : -howa  in  ? i rare  6.  Data 

was  recorded  on  tape.  Wat  t Lned  through  f ast-Fourier 

Detailed  descr  pt.i  ns  of  t •'.••id  e,.  . :ng  program  is  documented  by 


n 


analysis. 


SPECTRAL  DENSITY  ( M2  -SEC)  X 10 


i 


Figure  4 Comparisons  of  Numerical  Results  and  Analytical  Results  of 
Spectrum  Tranuforiii-.iti.cn  Over  Parallel  Hot  tom  Contour 
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Borchardt  (1974)  and  Hsu  (1975).  Samples  of  comparisons  are  shown  in  Figure  7. 
Since  there  is  no  simultaneous  measurement  of  deep  water  wave  conditions,  the 
deep  water  wave  spectrum  of  Pierson-Moskowitz  type  is  again  adopted  as  input 
condition  in  the  numerical  computation.  A north-easterly  wind  of  8.3  m/sec  is 
assumed  to  prevail.  Since  the  recorded  wind  speed  only  indicated  a range  of 
8-10  n/sec,  the  selection  of  8.3  m/sec  in  the  computation  is  merely  an  estima- 
tion, Varying  this  input  wind  speed  will  affect  basically  the  magnitude  of  the 
spectral  density  but  hardly  the  shape  of  the  function. 

Finally,  the  spectrum  transformation  over  a rythmic  nearshore  hydrograph 
such  as  the  type  described  by  Sonu  (1972)  is  demonstrated.  Assuming  that  the 
bathymetry  undulates  alongshore  with  diminishing  amplitude  toward  offshore,  as 
shown  in  Figure  8,  such  that  the  water  depth,  h,  in  the  x-y  plane,  can  be 
expressed  as  the  superposition  of  a uniformly  sloped  beach  and  a sinusoidal 
variation : 

h = -m(x-Ax)  + A (1  - *y)  sin  5-  (1  - -r*")  + h (20) 

o S 2 Ay  o 

for  x S and 
h = -m(x-Ax)  + hQ 
for  x > S 

where  m = beach  slope 

A = amplitude  of  shoreline  undulation 
o 

h = elevation  at  origin  of  coordinate  system 

o 

S = extent  in  the  off  hore  direction  where  oottom  undulation  begins, 
beyond  which  the  contour  is  parallel. 

Ax,  Ay  = grid  sizes  along  the  x and  y directions  respectively. 
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figure  8 Hypothetical  Rhythmic  Topogi  -phy 


Li 

lAi/~ 


¥ 

f 


-Convex  (A) 


Node  (B) 
-Concave  (C) 


B C 


Wave 


U - S.3  m/sec. 
0 = 0 


0.08  0.16  0.24  0.32  0.40  0.48 

FREQUENCY,  hz 


figure  9 Illustrations  of  Spatial  Variations  of  Wave  Spectrun 
Over  Phythriii  Topography 
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CHAPTER  VI  CONCLUSIONS 


A numerical  model  on.  wave  spectrum  transformation  over  irregular  bottom 
topography  has  been  developed.  This  model  proves  to  be  quite  effective  in  predicting 
wave  spectra  at  predetermined  locations  in  shoaling  water  outside  the  surf  zone 

provided  the  deep  water  wave  condition  is  specified.  The  numerical  scheme  has 
the  advantages  of  being  simple  and  stable,  in  particular,  for  cases  of  relatively 
mild  bottom  slopes.  The  model  is  linear  in  that  only  the  effects  of  shoaling 
and  refraction  are  considered.  Nonlinear  interactions  of  energy  transfer  from  one 
frequency  band  to  the  other  and  energy  dissipations  due  to  bottom  friction  and  tur- 
bulence are  neglected.  Thus,  its  application  to  water  of  very  shallow  depth  or 
near  the  surf  line  is  doubtful.  The  accuracy  of  results  for  cases  where  the 
bottom  topographies  vary  rapidly  at  the  side  boundaries  is  also  untested. 

The  model  has  been  verified  against  analytical  solution  for  the  case,  of 
parallel  bottom  contours.  As  to  the  comparison  between  the  predicted  result  and 
the  field  data,  the  agreement  is  considered  good  except  at  the  low  frequency  end. 
This  discrepancy  is  believed  to  be  mainly  due  to  the  assumed  deep  water  input  wave 
condition  rather  than  the  nonlinear  effects  mentioned  above.  Simultaneous  measure- 
ments of  deep  and  shallow  wave  conditions,  preferably  directional  spectrum,  should 
be  conducted  to  establish  the  field  applicability  of  the  model. 


It  should  be  noted  here  that  since  the  present  numerical  scheme  deals  with 
wave  height  rather  than  spectral  energy  density  as  the  basic  variables  in  the 
governing  equations,  the  incorporation  of  frictional  effects  and  turbulence 
energy  dissipation  is  a rather  simple  matter  if  such  dissipations  can  be 
empirically  related  to  the  diminishing  of  wave  height.  Furthermore,  the  inter- 
actions between  waves  and  currents  could  also  be  included  as  demonstrated  by  Nod 
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C A PROGRAM  TO  COMPUTE  SPECTRUM  TRANSFORMATION  OVER  IRREGULAR 

C BOTTOM  TOPOGRAPHY 

C GS  IS  ARRAY  OF  GRAIN  SIZE 

C OW  IS  ARRAY  OF  DEPTH  OF  WATER,  XBAR  IS  ARRAY  OF  SEA  LEVEL  INTERCEP 

C MI  = PAX  VALUE  FOR  I SUBSCRIPT,  NOT  TO  EXCEED  IDO 

C MJ  = MAX  VALUE  FOR  J SUBSCRIPT,  NOT  TO  EXCEED  1 HO 

C TITL  = IDENTIFYING  TITLE  FOR  EACH  SET.  (12A6). 

C NORAYS  = NO.  RAYS  IN  EACH  SET.  (15). 

C T = WAVE  PERIOD,  SECS,  (F8.2). 

C H = DEEP  WAIER  WAVE  HEIGHT,  (F6.2). 

C X , Y = STARTING  COORDINATES.  (F7.2). 

C A = INITIAL  DIRECTION  OF  RAY.  (F7.2). 

C 

COMMON  CU(5O,2(!),U(5O,2G),V<50,2O),Z<5O,2O),SI(5C,2C),CO(5O,2O>, 

1h( 50, 2C) ,CG(5u, 20, 5(50,20) ,H BREAK (50, 20), IB (50, 2U) ,DDDX (50,20) 

2,  D DD  Y ( C , c 0 ) ,C  (50,20) 

COMMON/CON/G,PI,PI2,RAD,EPS,DX,DY,DX2,DY2,T,SIGMA,M,N,IT 
COMMON  / S FRE  SS/  SI  GXY  (50, 20) 

D 1 MENS  I CN  QSX(1CC,1 GO ),QSY(1 00,100) 

DIMENSION  GS( 50,20) 

DIMENSION  TITL  (12),EMT(2) 

DIMENSION  SAVE  (50,20 
COMMON/DTI/DTI S(5C) 

CO  MMO  N / SPC / N 1 ( 20 >, S IG ( 2 C >, E 0 ( 2 0), E 1 ( 5 0,2 0) 

COMMON/CELW/DDU 
INTEGER  OPT 

97  FORMAT  (//’  THE  ORIGINAL  MATRIX  OF  WATER  DEPTH  IN  METER') 

908  F ORM AT  ( //, 1 X , ' T I D AL  LEV  EL=NN*  * , F6 . 2 ) 

98  FORMAT  (//,1X,‘  CX=*/E5.2,3X,  ,CY=',F5.2) 

99  FORMAT  (/////,  IX,  'N  1=  ’,  15, 3X,  'NJ= 15) 

100  FORMAT  (2  15)  ; 

101  FORMAT (1CX,7F1U.1/(SF10.1)) 

102  FORMAT  (3F5.2) 

103  F0RMATC15F  8.3) 

C READ  BASIC  DATA 

A= 1 3 5 . C 
NT  =5  . 

C NT-NO.  OF  DIVISION  OF  ONE  TIDAL  CYCLE,  NT  MUST  BE  ODD 

C NNT  = NO  CF  T ID  AL  CYCLES 

C NTS=  NO-.  OF  TIME  STEPS 

T P ER  = 1 2 .25 

N N T = 2 j 

DT=TPER  / FLOAT ( NT-1) 

N T S = N I * NN  T 
G=9.8C621 
P I =3 . 1 A 1 592  7 

P 1 2 = P I *2.0  j 

RAD= 180. C/PI 

READ  (5,100)  N I , N J 

C NI  IS  NUMBER  OF  GRID  POINTS  PERPENDICULAR  TO  THE  SHORE 

C NJ  IS  THE  NUMBER  OF  GRID  POINTS  PARALLEL  TO  THE  SHORE 

M=  NI 
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N = NJ 

I F (N  .6T.50.OR.NJ  .GT.20)  GO  TO  13 
READ  < 5 , 1C2) DX,  r'f  ,A3 

0 Xc=DX  *2.0 
DY  2=0  V*  c .0 
EPS=0  L C C 5 
NBOUN  0-5 
J J=N+2-NE0UND 
N P 1 - N ♦ 1 

READ  ( 5,  110)  ( < D V,  < I , J ),J=NBQUND,JJ),I=1,NI  ) 

FORMAT  ( 7 F 5 .1 > 

WRITE  (6,99).  N I , N J 
WRITE  (6,98)  D X , 0 Y 
WRITE  (6,908)  AC 
WRITE  (6,97) 

WRITE  (6,110)  ((  DW  ( I,J  ) , J=NB00N0,  J J ) ,1  = 1 ,M  ) 

DO  73  1-1, M 
DO  73  J=NRGUND,JJ 
DW(I,J)=DW(I,J)  + Au 
DO  15  1 1 = 1 , N I 
00  151  J=1,NR0FND-1 
DW  < I , J ) = DW (I ,N FOUND ) 

DO  1 5 G I = 1 , N 1 
DO  150  J = 1 2 , N J +1 

DW (I ,J) -CW( I ,J  J) 

WPITF  (6,108) 

WRITE  (6,103)  ( ( C U ( J,J),J  = 1,NPl),I=1,P) 

FORMAT  ( / / , 1 X , 'THE  AUGMENTED  MATRIX  OF  WATER  DEPTH  IN  METER') 

C NBOUND  JS  THE  NO.  CF  PARALLEL  PROFILES  POUNDING  THE  STUDY  AREA  ON  EACH 
THET  A 0 = A 

THETAO  = TI!ETAO/RAD 
THETA0=F12-THETA0 
T I M E - 0 . 

DO  13  K = 1 , 1 0 

SI G(K ) -L  .1*Pl * FLOAT  ( K) 

S 1 GM  A = S I C ( K) 

T = 2 . * P I / S I G ( K ) 

S I 4 = $ I G ( K ) * * A 
SI5-SI4*SIG(K) 

GL'=8 .3 

CO  NS  = CJ  ,7  4*  ( ( 9 . F/GU)  **4  ) 

AW=0.7779*EXPt-CONS/SI4)/SI5 
60 (K ) =0 . 1945 /CONS 

D R W~  F 9 ( K ) / ( A » ) 

HO -2  .<<28* SON T ( LO  (K) ) 

DO  14  1 1 T - 1 , 1 

CALL  RE F RAC ( THE T AO, HO, I T , N a 0 UND , J J , K ) 

106  CALL  LFUEL(NI,NJ,TIME) 

141  CONTINUE 
13  C 0 N T i N I , E 
STOP 

E N D 

... — 

min  i i i NMl 


1 1 c 


73 

151 

15G 


108 


p-r 

I 


SUBROUTINE  REFEAC(THETAO,HH,ITER,NBOUND,JJ,K) 

COMMON  C (5u,20>,U(50,2G),V<50,2U>,Z(50,20),SI(5f,2C),CO<50,20), 
lH<50,2C>,CG<5u,2C>,S<5U,2u>,HRREAK(5o,2U>,IB(50,,20)„DDDX(50,20> 
2,DDDY ( 5C  ,20)  ,C  (50, 20 

C0MM0N/CCN/G,PI,PI2,RAD,EPS,DX,DY,DX2,DY2,T,SIGMA,M,N 
COIMMON/CTI/DTI  S(50) 

COMMON/SFC/Hl(2G),SlG(2C),EO(2C>,El(50,20) 

COMMON/DELU/DDW 

K1=M-1 
Nl  =N+  1 
UP 1 = N + 1 
NBM1=NBCUND-1 
N2=N+  2 
CALL  DEPTH 

10  CALL  SNELL(THETAO,HH,NBOUND,JJ) 

CALL  ANCLE (20) 

CALL  H E IOHT ( 10 ,1  TER  ,K ) 

RETURN 

END 


SUBROUTINE  DEPTH 

COMMON  c (5Q,2G),U(50,2C),V(50,20),Z(50,20),SI(5P,2C),CO(50,20), 
lH(50,2C),CG(50,2G),w(50,20),HBRLAK(5C,2C),IB(50,2G),DDDX(50,20) 
2, ODD Y ( 5 C ,20 ) 

COMMON/CCN/G,PI,PI2,RAD,FPS,DX,DY,DX2,DY2,T,SIGMA,M,N,IT 
MM  1 = M- 1 

156  DO  157  1=2, MM1 
DO  157  J=2,N+1 
A = D ( I ♦ 1 , J ) 

IF  (A .LT  .f  .0)  A = 0 .0 
P = D ( 1-1, J) 

IF(B.LT.C.O)  B =0 .0 
C=  D( I , J * 1) 

IF(C.LT.C.O)  C = 0 . C 
E = D ( I,  J - 1 ) 

IF(E.LT.C.O)  E =0.C 
D D D X ( I,J)  = (A-B)/(2.0*DX) 

157  ODDY ( I ,J)  = (C-E )/ (2.G*DY  ) 

DO  mo  J = 1 , N + 1 

100  DDDX ( 1 , J)  = DD  DX  (2  , J) 

DO  101  1=1 ,M 

DDDY(  I,1)=DDDY  (1,2) 

101  DDDYl I ,N+2)=DDCY( J,N+1) 

158  RETURN 
END 
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SUBROUTINE  SNELL  ( TH  E T A 0 ,HH  , NB  0 UND  , J J ) 

COMMON  D <50,2U),U(50,2C),V<5C,2a),Z(50,20),SI(5C,20),CO<50,20>, 
1H(5U,2C),C6(5C,2C),S(SO,20),HPREAK(50,20),IB(50,20),DOOX(50,20) 

2 , D DO  Y ( 5C ,20 ) 

COMMON/CON/G,Pl,PI2,RAD,EPS,DX„DY,DX2,DY2,T,SIGMA,K,N 
DO  6 00  J = 1 , N ♦ 2 
30  DO  600  1 1=1  / 

I = M-  I I ♦ 1 

IF (D( I,J).GT .G.L)  GO  TO  33 
A N G=  P I 
WVHT  = 0 ;C 
SS=0 .0 
CC=-1  .C 
GO  TO  A3 
33  CONTINUE 

CALL  WVMM  ( D < I ,J  ) ,U  .Q,O.C,0 .0,G.0,RK  ,SIGPA) 

A A =R  K * D ( I , J ) 

ANG=ASIN(SIN(THETAO)*TANH(AA>) 

ANG=P1 -ANG 
ARG=2.0*AA 

SHOAL  = SQRT(1.G/(TAUH(AA)*(1.U  + ARG/SINH(ARG)))) 

RE  F = SQRT (COS (THET  AO) /COST ANG) ) 

WVHT=HH*SHOAL*FE  F 
H H =0 . 1 2 * T A N H ( A A ) * P I 2 / R K 
I F (WVHT-.GT  . H B ) WVHT  = H3 
I N = 1 

I F (WVHT-  .GT  .HP)  I N = 0 
SS  =S I N ( ANG ) 

CC=COS  (ANG) 

43  CONTI N U E 
U(I,J)-C.Q 
v ( i , j ) = r .o 

H ( I,J ) = * V H T 
IB (I , J ) = IN 
/ ( I , J ) = A N C 
SI ( l , J ) = SS 
CO (I , J) =CC 
60 U CONTINUE 
700  RF  TURN 
END 
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SUBROUTINE  GROLP(I,J,DCGDX,DCGDY,FF,RK) 

COMMON  C (50, 20)  ,U  C 50, 2G),V(50, 20 ), Z < 50,20) , 51(5  0,2  0,  CO  (50, 20), 
lH(50,2C),CG(5Li,2l),S(50,2Q),HBREAK(5G,2U),lB(5U,20),DDDX(50,20) 
2,DDDY(5C,20),C2(5C,20) 

COMMON/CON/G,PI,P12,RAD,EPS,DX,DY,DX2,DY2,T,SIGMA,M,N 
DIMENSION  DUD X(  50, 20)  „DUDY  < 5 0,20) ,DVDX (50,20 >, DVD Y (50,20) , 

3 DTDY(50,23),DKDX(5Ci,2G),DKDY(50,20),E(50,20),DTDX(50,20) 

JJ=J 

DE  P=D  ( I , JJ  ) 

COSI  = CO  ( 1 , J > 

S I NI  = S 1 ( l,j ) 

DCGDX-O.C 
DCGDY=O.C 
F F =0 . 0 

CALL  WVM'K(DEP,COSI,SINI,U(I,J),V(I,J),RK,A) 

T A =T  A N H (RK*DFP  ) 

HBREAK (I,J)=0.12*PI2*TA/RK 
C0SH1 =CCSH(RK*  CEP) 

SE  CHSQ=  1 .0/ (CO  SH 1**  2) 

ARG2  = 2 .'C*HK*DEF 

S 1 NH  2 = S 1 NH  ( A KG  2 ) 

COSH2=CCSH(ARG2) 

SINHSQ=SINH2** 2 

E(I,J>=U(I,J>*COSI*V(I,J>*SINI+0.5*A*(1„C*ARGZ/SINH2)/RK 
E E =E ( I,J) 

C - SO  R T (L*TA/RK  ) 

FF=U.5*  (1.0+ARG2/SINH2) 

CG  ( I , J)  = FF  *C 

DUDX(I,J)=(U(I+1,J)-U(I-1,J))/DX2 
DUDY ( I , J)  = (U( I ,J ♦ 1) -U( 1 ,J-1 ) ) /DY2 
DVDX(I,J)=(V(I+1,J)-V(I-1,J))/0X2 
DVDY(I,J)  = (V(I,JM)-V(I,J-1))/DY2 
DTDX  ( I , J ) = ( Z ( I ■*1,J)-Z(I  - 1 , J ) > / D X2 
DTDY(I,J>=(Z(I,J+1)-Z(I,J-1))/DY2 

DKDX(I,J)=RK*((U(I,J)*SINI-V(I,J)*C0SI>*DTDX(1,J)-(C0SI*DUDX(I,J)+ 
1SINI*DVDX(I,J)  )-A*DDDX(I,J  ) / S1NH2  > /EE 
DKDY(I,J)  = RK*(<U(1,J)*S1M-V(I,J)*C0SI)*DTDY(I,J)  -(COSI*DUDY(I,J) 
1 +SINI *CVDY( I, J))-A*DDDY(I,J  )/SJNN2)/EE 
P=C* ( S INH2-ARG  2*COSH2) / S I NH  SO 
DKDDX  = RIC*DDDX(  1 , J J ) +DEP*DKDX(  I , J ) 

DKDDY=RK*DDDY(  I,JJ)+DEP*DKDY<  I , J ) 

Q=0.5*G/(C*RK**2) 

C3 ( I ,J)  = C 

OCDX=Q*(RK*SECHSG*DKOOX-TA*OKDX(T,J)> 

DCDY  = Q*(RK*SF.  CHSG*DKDDY-TA*DKDY(I,J)) 

DCGDX=P*DKDDX+FF*DCPX 

DCGDY=P*CKDDY*fF*PCDY 

RETURN 

END 
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SUBROUTINE  LEV  EL  ( N 1 , N J , T 1 me  ) 

COMMON  CW(5D,20) , II  < 5Q,2C),V ( 50,20) , / (50,20) , S I ( 5P,2C) ,C0 (50,20), 
lF)<50,20,CG(5C,2:),S(5U,20),HBkEAK(5C,20,IB(50,20),DDDX(50,20> 
2, DOPY (51, 20) ,C  (5C,2C) 

C0MM0N/C0N/G,PI,PI2,RAD,EPS,DX,DY,DX2,DY2,T,SIGKA,M,N 
DIMENSION  XBAR  (20 

C LINEAR  INTERPOLATION  OF  SEA  LEVEL  POSITION 

111  FORMAT ( 'NO  SLA  LEVEL  INTERCEPT  ALONG  LINE’,  15) 

104  FORMAT  (7F 6.2) 

REAL  INCREX 


DO  2 J=1,NJ 
DO  3 1 = 1, NI 

1 F (DU  ( 1 , J) .GT . F .)  GO  TO  4 

3 CONTINUE 

WR  IT  E (6 , 11  1 ) J 
GO  TO  2 

4 INCREX  = -CX*DW(  T,J)/(DW(I-1,J)-DW(I,J)> 

XBAR (J ) = FLOAT  < I)*DX“DX/2, -INCREX 

2 CONTINUE 
RE  TURN 
END 


SUBROUTINE  NEVfHT  (I,J,1FLAG) 

COMMON  C (50,20) ,U(  50,20) ,V(  50,20) ,Z (50,20) ,S I < 50,20 , CO (50,20), 
1H  < 50, 20, C G(50, 2C ), S(50, 20), HBREAK(50, 20, IB(  50, 20), DDDX  (50, 20) 
2,DDDY ( 5 L ,20 ) 

COMMON/CON/G,P1,PI2,RAD,LPS,DX,DY,DX2,DY2,T,S1GMA,M,N 
C 0 mm on/xpp/x break  (20),SAVE(  50,2  0 
SAVE (M  , J)=H («, J) 

IB (I , J)  = 1 

Ml = N ♦ 1 
N2  =N  ♦ 2 

' CC1=(V(I,J)+CG(1,J)*SI(I,J))/DY 

CC2=(U(I,J)*CG(I,J)*C0(I,J))/DX 

HN  EW=  <CCl*H(I,J-1)-CC2*H(IM,J))/(CCl-CC2-S(I,J)/2.0) 

I F ( H N E k . LT.  C.O)  HN Erf  = 0.0 
SAVE ( 1 ,J)=HNEW 

I F (KNEW  .LE  .HPREAKd  ,J  ) ) GO  TO  850 
3 HNEW  = HF'REAK  ( 1 , J) 

IB (I ,J) =L 
850  CONTINUE 

860  I F (ABS  (HNEW-H ( I ,J  ))  .GT  . (EPS *ABS (HNEW) ) ) 1 F L A G = 0 
H( I,J)=HNEW 
899  RETURN 

END  ______ 

I 

I 
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SUBROUT  IKE  HEIGHT  (ITVAX,KK,K) 

COMMON  C (50,2C>,U(50,2U),V(50,20),Z(50,2G>,SI(5C,2C>,CO(50,20>, 
lH(50,2C),CG<50,2C),S(50,20),HBRfcAK(50,20),lB(50,20),DDDX(50,20) 
2,D0DY ( 5C#20) ,C (50,20) 

COMMON/CON/G,PI,PI2,RAD,EPS,DX,DY,DX2,DY2,T,SIGMA,M,N 

COMMON/STRESS/S1GXY(50,21j> 

COMMON  / X (•  R / X U R EAK  (20  ,SAV£  ( 5 G , 2 C) 

COMMON  / C Cl  J / CC  (50,20) 

COMMON/SPC/Hl(20),SIG(2C),tO(2C>,E1(50,20) 

COMMON / DELW/CDU 
DIMENSION  AA  W ( 50, 20 ) 

DIMENSION  AZ(5C,20) 

COMM ON / DTI/DTI  S(50) 

M 1 =M- 1 

N 1 =N ♦ 1 

C COMEPUTE  VALUES  OF  S(I,J> 

DO  500  1=2, Ml 
DO  5 U 0 J =2,N  1 

I F (D( J , J).LE .0.0)  GO  TO  500 

CALL  GROL'PC  I ,J  ,DCGDX,DCGDY,FF,RtO 

DUDX=(U(I+1,J)-U(I-1,J))/DX2 

DUDY=(U(I,J+1)-U(I,J-1))/DY2 

DVDX=(V(1+1,J)-V(I~1,J))/DX2 

D VDY=(V(I,J  + 1)  -V(I,J“1  ) ) / D Y 2 

DTDX-CZ  ( 1 + 1 , J ) ~Z  ( 1-1, J)  ) / D X 2 

DTDY=(Z(I,J+1)-Z(I,J-1))/DY2 

S S 2 = S I ( I , J ) * *2 

CC  2 = C0  ( I ,J ) **2 

S I GX  X = ( 2 .0*  F F-C'.5)*CC2  + (FF-0.5)*SS2 
SI  GY Y = ( 2.G*F F-C .5>*SS2+ (FF-0.5)*CC2 
TAUXY=FF*SI(I,J)*CO(I,J) 

SIGXY  ( I , J ) = T AUXY 

S<1,J)=CG(I,J)*(SI(I,J)*DTDX-C0(I,J)*DTDY)-(DUDX+DVDY)-(C0(I,J)*DC 
1GDX+SI(I,J)*DCGDY)-(SIGXX*DUDX+TAUXY*DUDY+TAUXY*DVDX+SIGYY*DVDY) 
500  CONTINUE 
N2  =N  + 2 
M2  =M-2 

DO  510  11=1, M2 
I =M-  I I 

DO  580  1T=1,1TMAX 
I F LAG  = 1 
DO  520  J=2,N1 

IF (D< 1,J).LE .0.0)  GO  TO  520 
CALL  NEWHT ( I ,J  ,1  FLAG) 

520  CONTINUE 

I F (I  FLAG  .EQ  .1)  GO  TO  570 
580  CONTINUE 

WRITE  (6,540)  I , I T 

540  FORMAT(*  RELAXATION  FOR  THE  WAVE  HEIGHT  FAILED  TO  CONVERGE*, 

TON  R0E,,IS,*AFTrR,,I6,'ITERATI0NS*) 

WRITE (6,541) (H(I,J),J=1,N2  ) 

541  FORMAT  (*  LAST  VALUES  Of  H A R E * /, ( 1 C F 1 3 . 5)  ) 

RETURN 

570  WR ITE ( 6 ,542)  I , I T 

542  FORMAT  ( 1 CX  , * RELAXATION  FOR  WAVE  HEIGFIT  CONVERGED  ON  R OW  * , I 6 , 3 X , * A 
IF  T ER  * , 1 6 ,* 1 TER  AT  I ONS* ) 

510  CONTINUE 
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WRITE  (6,6052) 

6052  FORMAT  ( / / , 1 X , ‘WAVE  HEIGHT  MATRIX'/) 
WRITE(6,652)((Fi<I,J),J  = 1,Nl),I  = 1,M) 

WRITE  (6,6053) 

6053  FORMAT  <//,  1 X , ’BREAKING  POSITION  ("0" MEANS 
WRITE(6,653)((IET(I,J),J  = 1,N1),I  = 1,M) 

653  FORMAT  (1518) 

652  FORMAT  ( 15F8.7) 

00  57  5 J = 1 , N 1 
A = 0 . 


DO  560  11=1, M2 


WAVE 


BREAKING  )'  ) 


I = M-  I 1 

IF  (A-.GT.O.)  C-0  TO  575 

o x x = n . 

DDXX=HEKEAK(l4l,J)-HBPEAK(I,J) 

IF  (SAVE  ( I,J)-HEREAK  ( I , J ) >560,560,554 
554  DDX=((SAVE(I,J)-HhREAK(I,J))/(SAVE(I,J)-SAVE(I-*1,J)*DDXX))*DX 


557  XBREAK(J)=FL0AT(I)*DX-DX/2.+DDX4DXX 
A=XBRE  A K ( J ) 

60  TO  5 6 C 
555  DDX=n. 

XBREAK(J)=FL0AT(I)*DX-DX/2.+DDX+DXX 
A = X BP  E AK ( J ) 

560  CONTINUE 
575  CONTINUF 
NBCUN  D = 5 
NB UNBOUND 
J J =1  1 

DO  702  I = 1 , M 

DO  702  J =NB , J J 

El  (I , J > =H< I , J) **2/8  . 

702  AAW( I , J ) = E 1 (I , J)  / (DDW) 

WRITE  (6,112)  SIGMA 

112  FORMAT  (//, IX,  'OFFSHORE  ENERGY  DENSITY  AT  S I GMA = ' , F 8 . 3 ) 
WRITE  (6,113)  EO(K)  ,AW 

113  FORMAT  (10X,'E,C(K)=',F1C.3,10X,'AW=',F10.3) 

WRITE (6,704) SIGMA 

704  FORM  AT ( / /,1 X , ' NE ARSHOR E ENERGY  DENSITY  AT  S I GM A = ' , F 8 . 3 ) 
WRITE (6,651) (( AAW(I , J ) ,J  = NB,JJ  ),I=1,M) 

651  FORMAT  ( /F10. 4) 


6051 


WRITE  (6,6051)  SIGMA 
FORMAT  (//,1X,*fcND  C F 
1 • ♦ 

R E TOR  N 
END 


CALCULATION  AT 


S I c,  M A = 
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SUBROUTINE  ANGLE(ITMAX) 

COMMON  D (5i>,20),U(50,2C),V(50,20)  ,Z  (50,20)  ,SI(  50,20  , CO  (50,20), 
lH(5U,2C>,CG<50,2C),S(5(j,2a),HBREAK(50,20),IB(50,20),DDDX(50,20) 

2,D  DD  Y < 5C  ,20) 

COMMON/CON/G,PI,P12,RAD,EPS,DX,DY,DX2,DY2,T,SIGMA,K,N 
Nl =N+ 1 
N 2=N  + 2 
Ml =M-1 
m2  =m-2 

DO  200  IT-1, UMAX 
1 F LA  G = 1 
DO  210  J=2,N1 
DO  211  1 1=2, Ml 

I =M- I 1 ♦ 1 

I F (D  ( I , J ) .LE  .O.C-.OR  .D  ( r-1,J  ) .LE  .0.0)  GO  TO  211 
CALL  N E V.  ANG  ( 1 , J , I FL  AG,  2 ) 

211  CONTINUE 
2 1 0 CONT I NUt 

I F (I FL AG  .EQ  .1)  6 0 TO  25  C 
200  CONTINUE 

WRITER,  220)  UMAX 

220  FORMAT ( * RELAXATION  FOR  THETA  FAILED  A F T E R * , I 6, 3 X , ' I T E R A T I 0 N S * ) 

250  WRITE(6,251)  I, IT 

251  FORMAT ( 1 LX , * RELAXATION  FOR  THETA  CONVERGED  ON  R OU • , I 6 , 3 X , * A 
IFTER',16,*  ITERATIONS*) 

WRITE  (6,2501) 

2501  FORMAT  (//, IX,  ’WAVE  APPROACHING  ANGLE  IN  RADIAN*) 

WRITE (6,252) ( ( 2 ( I ,J),J=1,Nl),I=1,M) 

252  FORMAT  (15F8.3) 

RETURN 

END 
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SUBROUTINE  NEW ANG(I ,J, IFLAG,I M ) 

COMMON  C <50,2C),U(50,2G),V  <50,20>,Z (50,20) ,SI (50,20 , CO (5 0,20), 
lH(5o,2C),CG<50,2C>,5<50,20),HBREAK(S0,2G),IB(5C,2U),DDDX(50,20) 

2,  D DD  Y ( SC, 20) 

CO  MM  ON/C  ON  /G,P1,PI2,  RAD,  EPS,  DX,DY,DX2,DY2,T,  SIGMA,  M ,N 
C STATEMENT  FINCTIONS 

DUDX(I,J)=<U<1*1,J)-U(I-1,J))/DX2 
OUCY ( 1 ,J)  = (U(I  ,J  + 1) -U(I  , J - 1 ) ) / PY  2 
DVDX(I,J)=(V(I+1,J)-V(I-1,J))/DX2 
DVDY(I,J)=(V(I,J+1)-V(I,J-1))/DY2 

F ( 1, J ) = 11 (I , J )* COS  I ♦ V( I , J )*S IN  I ♦ 0.5*A*(1.J  + ARG2/SINH2)/RK 
DKDY(I,J)=(-(CCS1*PUDY(I,J>  ♦ S1NI *DV DY ( I, J ) ) -A* DD DY( I ,J  )/SINH2 
1/F  F 

DKDX(I,J)=(-(CCS1*DUDX(1,J)+SIN1*DVDX(I,J))-A*DDDX(I,J  > / S INH2) / 
IF 

FAC( I , J )=U( I ,J  )*SINI-V( I , J ) * C 0 S I 
COSI -CC ( I, J ) 

S I N 1 = S 1 < I , J > 

J J =J 

CALL  WVNLMC  D ( I ,JJ),C0SI,$1NI,U(I,J),V(I,J),RK,A> 

ARG2  = 2 .C*RK*P(  I,JJ) 

SINH2  = S JNH( ARG2) 

FF=F(I,J) 

I F (F  F .GT-.O.C)  CO  TO  450 

WRITE  ( 6 , 4 S 1 ) I , J , D ( I,  J J ) ,COS  I ,SIN1,U<1,J),V(I,J)  , R K ,A 
451  FORM  AT  ( 1CX  , * f F IS  NEGAT  IVE  — OUTPUT  I , J , D , CO  S I ,S  I N I ,U,  V ,P  K , A * / , 

12  I 5, 7 El  3 .4) 

RETURN 

4 SC  F AC1  = F AC (I , J ) 

DENl  = (SJM-COSI*FACI/FF>/DY 
DEN2=(CCSI+SINI*FAC1/FF  > /DX 
PE  N=0EM-DEN2 

ZNEW=<CCSI*DKDY(I,J)-SINI*DKf>X(I,J)+Z(I,J-1)*DENWU-H,J)*DEN2)/ 

1EN 

IF  (AF<S(ZNFW-Z(  I , J ) ) .GT.(EPS*ABS(ZNEW)>)  I F L A G = 0 

l ( I, J ) =ZNFW 

CO ( I , J) = COS ( 7 ( I, J) ) 

SI (I,J>=S1N(Z(I,J>> 

RETURN 

END 


I 

l 
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SUBROUTINE  WVNUM(D,C0SI,SINI,U,V,RK,A) 

SUBROUTINE  COMPUTES  THE  WAVE  NUMBER  INCLUDING  W A V E-  C OR  R EN  T INTERATION 
C0MM0N/C0N/G,PI,PI2,RAD,EPS,DX,DY,DX2,DY2,T,SIGMA,M,N 


EP  SK  = 2 .CC05 
R K =P  1 2 / ( T*  SO  RT  ( G * D ) ) 

DO  10U  1=1, 50 

A=SI6MA-ll*RK*CCSI-V*RK*SINl 
A2  =A * * 2 
A H G = k K * t 
F 1 =EX  P (ARC,) 

F 2 =1  .0/  FI 

S£CH  = 2'.C/(FlfF2) 

SECH2=SECH*SECH 
T T =T  A N H (ARG  ) 

FK=G*RK*TT-A2 

FFK=G*(ARG*SECH2yTT)42.C*(U*C0SI+V*SINI)*A 
RKNEW=RK-FK / F F K 

1 F (ABS ( KKNEW-R K)  .LE  .( ABS (EPSK*R<NEW)  ) > GO  TO  110 
RK=RKNEW 

100  CONTINUE 

W R IT  E ( 6 , 1G1  ) I,RK,T,D,U,V 

101  FORMAT!'  ITERATION  FOR  Y FAILED  TO  CONVERGE  AFTER*  ,I6,3X,5FlO.O 
RETURN 

1 1 C R K =R  K N F W 

A-$IGMA-U*RK*COSI_V*RK*SINI 
IF  CRK  .GT  .0.0)  CO  TO  120 

WRITE (6,  130)  D,COS),SINI,U,V,RK,A 

130  FORMAT ( 1CX,*RK  IS  N EGA T I VE--OUTPUT  D , C 0 S I , S I N I , U , V , RK , A * / , 7 El  2 .4 ) 
120  RETURN 
END 


WAVE  APPROACHING  ANGLE  IN  RADIAN 
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ILLUSTRATIVE  SAMPLE 


COMPUTATIONS 


riuc-D 


The  computer  program  has  also  been  converted  to  suit  the  ICL  1906 
machine  at  the  Technical  University  at  Braunschweig,  West  Germany.  Field  work 
has  been  conducted  at  the  Island  of  Sylt  in  an  attempt  to  verify  the  model. 

A set  of  sample  computations  are  provided  here.  Data  reduction  is  in  progress 
and  the  results  will  be  reported  separately  at  a later  date. 
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FIGURE  10  Location  of  Test  Site  (Island  of  Syl t/Germany , From  Dette,  1974) 
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Based  on  the  assumption  that  the  wave  energy  associated  with  a narrow 
frequency  band  stayed  within  the  band  upon  refraction,  a numerical  model 
based  on  finite  difference  method  has  been  developed.  This  model, 
utilizing  electronic  computer,  proves  to  be  quite  convenient  to  determine 
shallow  water  wave  spectra  at  designated  spatial  locations  of  irregular 
bottom  contours  provided  the  deepwater  wave  spectrum  is  given.  The  model 
has  been  verified  against  analytical  solution  for  case  of  two-dimensional 
parallel  bottom  contours  and  compared  in  good  agreement  with  field  measure- 
ment at  location  well  beyond  the  surf  zone. 
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